CAM Implementation of Mesoscale Eddy Diffusivity in Terms of Mixing Length Theory

MICHAEL BATES, ROSS TULLOCH, JOHN MARSHALL, AND RAFFAELE FERRARI

Steering Level status

1. L_eddy values ok
2. Rossby wave speed ok
3. Eady inverse timescale 
    a.has anomalous equatorial belt (f goes to 0 at equator)
    b.Richardson averages are too high causing sigma to drop too low
4. Zonal Eddy Phase Speed (c) 
    a. Too high 10S:10N
    b. Missing westward velocities in ACC which appear 
        in Hughes Data used in this parameterization
5. (U-c) Because of problem with c this term is also too high 10S:10N.
6. U_rms  (alpha*sigma*l_eddy)        
    a. Overall too weak  (due to problems with Eady timescale.)
    b. Missing alot of structure especially in the NH. Western Boundary currents missing
        some wbc structure in eady but multiplying by leddy removes that at higher latitudes
7. Zonal mean velocity looks ok compared to ECCO annual average
8. Suppression, Lmix, and K are wrong because of errors in the above terms

Todo

1. New calculation for sigma may correct sigma deficiencies and fix the U_rms problems.
    or fixing equator with max (/f/,sqrt(c_r*2*Beta)) and changing calculation for sqrt(R_i) would improve sigma
2. Need to work on Zonal phase speed c
3. Fill in rest of plots for K
4. Plot vertical distribution of K

In [1]:

%pylab inline
from IPython.display import display
from IPython.display import Image
from netCDF4 import Dataset
from netCDF4 import num2date
import numpy as np
import sys
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
from matplotlib.colors import LogNorm   
from matplotlib.ticker import LogFormatter
from mpl_toolkits.basemap import Basemap, shiftgrid, addcyclic
#for interpolation (you will have to install pyresample first)
matplotlib.rc('text', usetex=True)
matplotlib.rcParams['text.latex.preamble']=[r"\usepackage{amsmath}"]

import pyresample

numpy.version.version

Populating the interactive namespace from numpy and matplotlib

Out[1]:

'1.9.1'

In [2]:

def make_basemap(projection='cyl', resolution='c',lon_0=0, ax=None,llcrnrlon=25):
    m = Basemap(projection=projection, resolution=resolution,
                ax=ax,llcrnrlat=-90,urcrnrlat=90,
                llcrnrlon=llcrnrlon,urcrnrlon=llcrnrlon+360,fix_aspect=False)
    m.drawcoastlines()
    m.drawmapboundary(fill_color='white')
    m.fillcontinents(color='k',lake_color='white')
    parallels = np.arange(-60, 90, 30.)
    meridians = np.arange(-360, 360, 60.)
    m.drawparallels(parallels, labels=[1, 0, 0, 0])
    m.drawmeridians(meridians, labels=[0, 0, 1, 0])
    return m

def MapZoneContour(lon,lat,data,figsize=(14,9),title=None,levels=None,fmt=None,norm=None):
  
    map2zonal=6
    
    fig = plt.figure(figsize=figsize)  
    
    ax1 = plt.subplot2grid((1,map2zonal), (0,0), colspan=(map2zonal-1))
    
    m=make_basemap()
    
    #make longitudes monotonic for contouring add longitude period to end
    lon_mono = np.where(np.greater_equal(lon,min(lon[:,0])),lon-360,lon)+360.
    endlon=lon_mono.shape[1]-1
    lon_mono[:,endlon]=lon_mono[:,endlon]+360.

    # contour and labels
    if (levels==None):
        CS=m.contour(lon_mono,lat_targcyc,data,10, linewidth=.5, colors='k')
    else:
        CS=m.contour(lon_mono,lat_targcyc,data,levels=levels, linewidth=.5, colors='k')
        
    plt.clabel(CS, CS.levels,fmt=fmt)
    plt.suptitle(title,fontsize = 15,y=1.04)

    # zonal subplot
    ax2 = plt.subplot2grid((1,map2zonal), (0,(map2zonal-1)), sharey=ax1)
    
    # First find the zonal mean by averaging along the latitude circles
    data_zonal = data.mean(axis=1)

    plt.plot(data_zonal,lat[:,0])
    return fig

In [3]:

def MapContourf(lon,lat,data,
                figsize=(14,9),
                title=None,
                levels=None,
                fmt=None,
                extend='neither',
                contourlines=False,
                addzonal=False,
                invertcb=False,
                norm=None,
                yscale='Linear'):
#    from matplotlib.colors import LogNorm   
    map2zonal=6

    formatter = ticker.ScalarFormatter()
    formatter.set_scientific(True)
    formatter.set_powerlimits((0,0))

    fig = plt.figure(figsize=figsize)  

    #make longitudes monotonic for contouring add longitude period to end
    lon_mono = np.where(np.greater_equal(lon,min(lon[:,0])),lon-360,lon)+360.
    endlon=lon_mono.shape[1]-1
    lon_mono[:,endlon]=lon_mono[:,endlon]+360.
   
    if (fmt==None):
        fmt='%.3e'
    
    if (addzonal):
        ax1 = plt.subplot2grid((1,map2zonal), (0,0), colspan=(map2zonal-1))
    else:
        ax1 = fig.add_subplot(111)

        
    if (contourlines):
    # Create contour lines
        CS0 = ax1.contour(lon_mono,lat,data,levels=[0.],norm=norm,linewidths=0.5,colors='k',latlon=True)
        
    # y-labels on the right
    ax1.yaxis.tick_right()
      
    if (invertcb): 
        CS1 = ax1.contourf(lon_mono,lat,data,levels=levels,norm=norm,cmap='jet_r')         
        cbar = fig.colorbar(CS1, ax=ax1,format=fmt)
        cbar.ax.invert_yaxis()
        if (extend!=None):
            CS1.cmap.set_under('darkred')
            CS1.cmap.set_over('black')
    else:
        if (norm==None):
            CS1 = ax1.contourf(lon_mono,lat,data,levels=levels,extend=extend,norm=norm,cmap='jet') 
            cbar = fig.colorbar(CS1, ax=ax1, format=fmt)
        else:
            CS1 = ax1.contourf(lon_mono,lat,data,levels=levels,norm=norm,cmap='jet') 
            l_f = LogFormatter(10, labelOnlyBase=False)
            cbar = fig.colorbar(CS1, ax=ax1, ticks=levels, format=l_f)
        if (extend!=None):
            CS1.cmap.set_under('black')
            CS1.cmap.set_over('darkred')

    m=make_basemap()


    plt.suptitle(title,fontsize = 15,y=1.04)
    fmt = matplotlib.ticker.FormatStrFormatter("%.3f")
       
    if (addzonal):
        # zonal subplot
        ax2 = plt.subplot2grid((1,map2zonal), (0,(map2zonal-1)), sharey=ax1)
        formatter1 = ticker.ScalarFormatter()
        formatter1.set_scientific(True)
        formatter1.set_powerlimits((0,0))
        ax2.xaxis.set_major_formatter(formatter1)
        # First find the zonal mean by averaging along the latitude circles
        data_zonal = data.mean(axis=1)
        if (levels!=None):
            pylab.xlim([levels.min(),levels.max()])
        plt.plot(data_zonal,lat[:,0])
    return fig

In [34]:

basedir='/scratch/jet'

### get Variables from file

#nc = Dataset(basedir+'/steering/g.e11.GIAF.T62_gx1v6.steer.004.pop.h.nday1.0249-01-02.nc')
#nc = Dataset(basedir+'/steering/g.e11.GIAF.T62_gx1v6.steer.006.pop.h.nyear1.0249.nc')
nc = Dataset(basedir+'/steering/g.e11.GIAF.T62_gx1v6.steer.007.pop.h.nday1.yr.0249.nc')
#print(nc.variables)
sigma_avg = nc.variables['SIGMA_AVG'][0,:,:]
utgrd_max = nc.variables['UTGRD_MAX'][0,:,:,:]
vtgrd_max = nc.variables['VTGRD_MAX'][0,:,:,:]
stgrd_max = nc.variables['STGRD_MAX'][0,:,:,:]
ttgrd_max = nc.variables['TTGRD_MAX'][0,:,:,:]
###sigma_gamma = nc.variables['SIGMA_GAMMA'][0,:,:]
###maxftinv = nc.variables['MAXFTINV'][0,:,:]
fcor = nc.variables['FCOR'][0,:,:]
c_eddy = nc.variables['C_EDDY'][0,:,:]
c_rossby = nc.variables['C_ROSSBY'][0,:,:]
l_eddy = nc.variables['L_EDDY'][0,:,:]
l_mix = nc.variables['LMIX'][0,:,:]
l_rhines = nc.variables['L_RHINES'][0,:,:]
l_rossby = nc.variables['L_ROSSBY'][0,:,:]
l_rossbyeq = nc.variables['L_ROSSBYEQ'][0,:,:]
urms_avg = nc.variables['URMS_AVG'][0,:,:]
urmsavgsq = nc.variables['URMS_AVG_SQ'][0,:,:]
uminuscsq = nc.variables['U_MINUS_C_SQ'][0,:,:]
suppress = nc.variables['SUPPRESS'][0,:,:]
umean =  nc.variables['U_MEAN'][0,:,:]
umean10 =  nc.variables['U_MEAN'][0,:,:]
#umean50 =  nc.variables['U_MEAN_50'][0,:,:]
#umean100 =  nc.variables['U_MEAN_100'][0,:,:]
#umean200 =  nc.variables['U_MEAN_200'][0,:,:]
#umean500 =  nc.variables['U_MEAN_500'][0,:,:]
kappalat =  nc.variables['KAPPA_LATERAL'][0,:,:]
#tinv =  nc.variables['TINV'][0,:,:]
grate =  nc.variables['GRATE'][0,:,:]
btp =  nc.variables['BTP'][0,:,:]
hmxl2 =  nc.variables['HMXL_2'][0,:,:]
xmxl2 =  nc.variables['XMXL_2'][0,:,:]
lat_orig = nc.variables['TLAT'][:]
lon_orig = nc.variables['TLONG'][:]
###time = nc.variables['time']
dz = nc.variables['dz'][:]
dzw = nc.variables['dzw'][:]
nc.close()
nc_ecco = Dataset(basedir+'/steering/ecco_ann_zonal.nc')
ecco_u=nc_ecco.variables['U'][0,0,:,:]
lat_ecco = nc_ecco.variables['LATITUDE_T'][:]
lon_ecco = nc_ecco.variables['LONGITUDE_UN179_181'][:]
nc_ecco.close()

###time_conv = num2date(time, time.units)

# this is the lefthand longitude starting point for plotting
llcrnrlon=25
# proportion of plot to side zonal average ie. 6:1
map2zonal=6

### use time_conv for plots using time as an axis
###fig = figure()
###plot(time_conv, air[:,10,10])
###fig.autofmt_xdate()
###  or ###
###contourf(lat[:],time_conv,air[:,:,10])
###colorbar()

In [35]:

# setup for using pyresample to interpolate plot

# make sure lon is in range of -180:180 for pyregrid
lon_orig = (lon_orig + 180) % (360) - 180 

# open dataset for output grid get latout,lonout
fc = Dataset(basedir+'/steering/co2flux_fossil_1751-2006-monthly_0.9x1.25_c20100204.nc')
lon= fc.variables['lon'][:]
lat = fc.variables['lat'][:]
fc.close()

# get 2D versions of the lat and lon variables add longitude start point here!
lon_targ, lat_targ = np.meshgrid(lon+llcrnrlon, -1.*lat)

# add cyclic point to output lat,lon arrays for plot routines
lon_targcyc = np.concatenate((lon_targ, lon_targ[:, 0:1]), axis=1)
lat_targcyc = np.concatenate((lat_targ, lat_targ[:, 0:1]), axis=1)

# make sure lon is in range of -180:180 for pyregrid
lon_targcyc = (lon_targcyc + 180) % (360) - 180

#When all the necessary inputs are read in memory, that's where we start using pyresample's magic.
#Create a pyresample object holding the origin (POP) grid:

orig_def = pyresample.geometry.GridDefinition(lons=lon_orig, lats=lat_orig)

#Create another pyresample object for the target (curvilinear) grid:
targ_def = pyresample.geometry.GridDefinition(lons=lon_targcyc, lats=lat_targcyc)
print len(lon_targcyc[0,:])
print len(lat_targcyc[:,0])

289
192

In [36]:

omega=7.2921e-5
deg2rad=pi/180
f=2*omega*sin(lat*deg2rad)

Update the matplotlib configuration parameters:

mpl.rcParams.update({'font.size': 18, 'font.family': 'serif'})

Implementation Notes on Diffusivity with Steering Level Suppression

12 November 2014

Starting with

$$K=u_{rms}∗L_{mix}$$

the general form for diffusivity, $K$, is given by equation (6) of Bates et al. (2014). Namely,

$$L_{mix} = {\Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - c|^2 /u_{rms}^2 (z=0)} $$

Here,

$u_{rms}$ is the root-mean-square (rms) eddy velocity;
$L_{mix}$ is the mixing length;
$L_{eddy}$ is the eddy diameter (depth independent);
$u_{mean}$ is the mean zonal velocity (resolved);
$c$ is the zonal eddy phase speed (depth independent);
$\Gamma = 0.35$;
$b1 \sim 4$.

Note that

$L_{mix} = {\Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - c|^2 /u_{rms}^2 (z=0)}$

Eddy Length Scales:

Rossby deformation radius = $L_r = {c_r \over |f|}$

Equatorial Rossby deformation radius $= L_{req} = {\sqrt{c_r \over 2\beta}}$

Rhines scale $= L_{Rh} = {\sqrt{u_{rms} \over \beta}} \sim {\sigma_{vi} \over \beta}$

Where, $c_r$ is the first baroclinic wave speed computed following equation (2.2) of Chelton et al. (1998) with $m=1$;

$f$ is the Coriolis parameter;

$\beta$ is the latitudinal variation of the Coriolis parameter; and

$\sigma_{vi}$ is the Eady growth rate given by

$\sigma_{vi} = {f \over \sqrt{R_i}}$

with $R_i$ the vertically integrated (over the 100 – 2000 m depth range) Richardson number.

So, any one these length scales could be used as an eddy length scale. An alternative is

$L_{eddy} = min (L_r, L_{req}, L_{Rh}).$

Eddy Velocity:

$u_{rms} = alpha*\sigma_{vi}*L_r$

where $\sigma_{vi}$ is the Eady growth rate based on local Richardson number and $\alpha$ is a scaling constant.

Zonal phase speed:

$c = - \beta * L_r^2$

Bates, M., R. Tulloch, J. Marshall, and R. Ferrari, 2014: Rationalizing the spatial distribution

of mesoscale eddy diffusivity in terms of mixing length theory. J. Phys. Oceanogr., 44,

1523-1540, doi: 10.1175/JPO-D-13-0130.1.

Chelton, D. B., R. A. deSzoeke, M. G. Schlax, K. E. Naggar, and N. Siwertz, 1998: Geographical

variability of the first baroclinic Rossby radius of deformation. J. Phys. Oceanogr., 28,

433-460.

Tullock, R., J. Marshall, and K. S. Smith, 2009: Interpretation of the propagation of surface

altimetric observations in terms of planetary waves and geostropic turbulence. J.

Geophys. Res., 114, C02005, doi: 10.1029/2008JC005055.

Questions:

In the 2D implementation, $u_{rms}$ and $u_{mean}$ specifications: upper-ocean vertically or integrated or at $z = 0$?

A: $U_{rms}$ is not depth dependent; both $u_{rms}$ and $u_{mean}$ are for surface only
In the 2D implementation, vertical profile will be specified by $N2(z)$?

A: Yes, $N^2(z) \over N_{ref}(z)$ to be more precise
Local $R_i$ use imbedded in sigma in $u_{rms}$ calculation?

A: No, $u_{rms}$ is depth independent
alpha = ?

A: trial and error
Cancellation of $f$’s in $u_{rms}$ calculation?
Zonal phase speed equation correct? Both $\beta$ and $L_r$ will be positive, producing $c < 0$ always. This appears to be in contrast with Tullock et al. (2009).

A: What is plotted in Bates is (U-c)

An alternate derivation of $\sigma_{vi}$ is:

$ R_i = {N^2 \over { ( \frac {\partial u}{\partial z} )^2 + ( \frac {\partial v}{\partial z} )^2} }$

$ N^2 = {{-g \over \rho_0 }\frac {\partial \rho}{\partial z}} $

After hydrostatic and geostrophic approximations

$f \frac {\partial v}{\partial z} = {{-g \over \rho_0 }\frac {\partial \rho}{\partial x}}; \quad\quad f \frac {\partial u}{\partial z} = {{g \over \rho_0 }\frac {\partial \rho}{\partial y}} $

$\frac {\partial v}{\partial z} = {{-1\over f}{g \over \rho_0 }\frac {\partial \rho}{\partial x}}; \quad\quad \frac {\partial u}{\partial z} = {{1\over f}{g \over \rho_0 }\frac {\partial \rho}{\partial y}} $

$\therefore$

$ R_i = {f^2 N^2 \over \underbrace{{ {g^2 \over \rho_0^2 } ( \frac{\partial \rho}{\partial y})^2 + {g^2 \over \rho_0^2 } ( \frac{\partial \rho}{\partial z})^2} }_\text{m^4}}\quad\quad = \quad\quad {f^2N^2 \over m^4} $

$\sigma_{vi} = {f \over \sqrt{R_i}}\quad = \quad {f \over \sqrt{{f^2N^2 \over m^4}}}\quad = \quad {{\cancel f m^2} \over \cancel f N}$

$RX_1 = RX_{east} = \Delta\rho_x = \rho_{i+1,j} - \rho_{i,j}$

$RY_1 = RY_{north} = \Delta\rho_y = \rho_{i,j+1} - \rho_{i,j}$

$RZ_1 = RZ_{k+1} = \Delta\rho_z = \rho_{k} - \rho_{k+1}$

$\displaystyle{1 \over L_{R_i}} \displaystyle\int_{2000m}^{100m} \left\lbrace { {-d\over\rho_0}{\frac {\partial \rho} {\partial z}} \over { {g^2 \over \rho_0^2 } \left[( \frac{\partial \rho}{\partial y})^2 + ( \frac{\partial \rho}{\partial z})^2\right] } } \right\rbrace dz$

Note: missing $f^2$ which will be cancelled when forming $\sigma_{vi}$

$\quad\quad$ so $\cdots$ this is not $R_i$

Implementation notes

For Level K

Numerator : Top $= -grav * RZ_{SAVE}(\cdots k+1) * dzwr(k)$

Denominator :

$ \begin{align} work1 = p25 & * ( RX(..,i_{east},k)^2 \\ & + RX(..,i_{west},k)^2 \\ & + RX(..,i_{east},k+1)^2 \\ & + RX(..,i_{west},k+1)^2 ) / DXT(i,j)^2 \\ \end{align} $

$ \begin{align} work2 = p25 & * ( RY(..,j_{north},k)^2 \\ & + RY(..,j_{south},k)^2 \\ & + RY(..,j_{north},k+1)^2 \\ & + RY(..,j_{south},k+1)^2 \\ \end{align} $

$work3 = {\left( TOP \over (grav^2*(work1+work2))\right)}*dzw(k)$

Notes:

1)Need to be careful at top and bottom of ocean
2)Accurate dzw(k) for each (i,j) to form $L_{R_i}$
3) When constructing $sigma$ itself, use $RZ_{SAVE}$ with a minimum N value
4) use eps2

CESM First Baroclinic Wave Speed $c_r \cancel {= {1\over\pi} \int_{-2000}^{-100} N(z)dz}$ Now it is over full model depth

In [37]:

#Resample (aka re-project, re-grid) the NCEP data to target grid. First with nearest neighbour resampling...
c_rossby_nearest = pyresample.kd_tree.resample_nearest(orig_def, c_rossby, \
        targ_def, radius_of_influence=500000, fill_value=None)

fig1=MapZoneContour(lon_targcyc,lat_targcyc,c_rossby_nearest,figsize=(14,9),
                    title="CESM First Baroclinic Wave Speed (cm/s)",
                    fmt='%.1i' )
c_r_chelton=Image(filename=basedir+'/steering/chelton_sfig1.jpg')
display(c_r_chelton)

First Baroclinic Rossby Radius Chelton vs CESM

CESM Rossby deformation radius = $L_{eddy}$

$L_r = {c_r \over |f|}$

$L_{req} = {\sqrt{c_r \over 2\beta}}$

$\require{cancel} \cancel {L_{Rh} = {\sqrt{u_{rms} \over \beta}} \sim {\sigma \over \beta}}$

$L_{eddy} = min (L_r, L_{req}, \cancel{L_{Rh}}).$

In [38]:

#####Resample (aka re-project, re-grid) the NCEP data to target grid. First with nearest neighbour resampling...
####l_eddy_nearest = pyresample.kd_tree.resample_nearest(orig_def, l_eddy, \
####        targ_def, radius_of_influence=500000, fill_value=None)

fcor_nearest = pyresample.kd_tree.resample_nearest(orig_def, fcor, \
        targ_def, radius_of_influence=500000, fill_value=None)
absf=np.abs(fcor_nearest)
print np.max(absf),np.min(absf)
l_rossby_calc=c_rossby_nearest/absf
print l_rossby_calc.max(),l_rossby_calc.min()

btp_nearest = pyresample.kd_tree.resample_nearest(orig_def, btp, \
        targ_def, radius_of_influence=500000, fill_value=None)
l_rossbyeq_calc=np.sqrt(c_rossby_nearest/(2*btp_nearest))
l_eddy_calc=np.minimum(l_rossby_calc,l_rossbyeq_calc)

### convert regridded radius from cm to km
cm2km=1.e-5
print np.min(l_eddy_calc)
l_eddy_calc_km=l_eddy_calc*cm2km
clevs=[10,20,30,40,50,60,80,100,150,230]
fig2=MapZoneContour(lon_targcyc,lat_targcyc,l_eddy_calc_km,figsize=(14,9),
                    levels=clevs,
                    title="CESM First Baroclinic Rossby Radius (km) $L_{eddy} \sim min(l_r,l_{req})$",
                    fmt='%.1i' )
l_r_chelton=Image(filename=basedir+'/steering/chelton_fig2.jpg')
display(l_r_chelton)

0.0001436 1.01893e-06
3.44952e+08 7234.85
7234.85

$\sigma$ Eady growth rate

$\sigma$ is the Eady growth rate given by

$\sigma = {f \over \sqrt{R_i}}$

Calculated in the model as $\cancel {\sigma = {max(|f|,\sqrt{c_r2\beta}) \over \sqrt{R_i}}}$ Just using $\sigma = {f \over \sqrt{R_i}}$

with $R_i$ the vertically integrated (over the 100 – 2000 m depth range) Richardson number.

In [42]:

##### Plot max(sqrt(c_r*2*beta),|f|)
orig_ecco_def = pyresample.geometry.GridDefinition(lons=lon_orig, lats=lat_orig)

####grate_nearest = pyresample.kd_tree.resample_nearest(orig_def, grate, \
#####        targ_def, radius_of_influence=500000, fill_value=None)
####c_rossby_nearest = pyresample.kd_tree.resample_nearest(orig_def, c_rossby, \
####        targ_def, radius_of_influence=500000, fill_value=None)
####btp_nearest = pyresample.kd_tree.resample_nearest(orig_def, btp, \
####        targ_def, radius_of_influence=500000, fill_value=None)
####fcor_nearest = pyresample.kd_tree.resample_nearest(orig_def, fcor, \
####        targ_def, radius_of_influence=500000, fill_value=None)
####tinv_nearest = pyresample.kd_tree.resample_nearest(orig_def, tinv, \
####        targ_def, radius_of_influence=500000, fill_value=None)
####
####tinv_calc=np.sqrt(2*btp_nearest*c_rossby_nearest)


####maxf=np.maximum(tinv_calc,np.abs(fcor_nearest))
absf=np.abs(fcor_nearest)
fig3=MapContourf(lon_targcyc,lat_targcyc,absf,
                    figsize=(14,5),
#                    title=" maxf = $max(|f|,\sqrt{c_r*2*B})$")
                    title="$|f|$")
##### Plot sqrt(R_i)
grate_nearest = pyresample.kd_tree.resample_nearest(orig_def, grate, \
        targ_def, radius_of_influence=500000, fill_value=None)
print  grate_nearest.min(),grate_nearest.max()
#clevs=np.linspace(0,2200,23)
clevs=np.logspace(0,6,20)
norm=LogNorm()
sqrtgrate_nearest=np.sqrt(grate_nearest)
fig4=MapContourf(lon_targcyc,lat_targcyc,sqrtgrate_nearest,
                 levels=clevs,
                 norm=norm,
                 figsize=(14,5),
                 title="$\sqrt{R_i}$ (log scale)")


##### Plot max(sqrt(c_r*2*beta),|f|)/sqrt(R_i)

sigma_avg_calc_day=(absf/sqrtgrate_nearest)*86400.

clevs=np.logspace(-3,0,30)
norm=LogNorm()
fig1=MapContourf(lon_targcyc,lat_targcyc,sigma_avg_calc_day,
                levels=clevs,
                norm=norm,
                figsize=(14,5),
                extend='both',
                title="Eady inverse time scale calc $(days^{-1}) \sim {|f| \over \sqrt{R_i}}$ (log scale)")

clevs=arange(.005,.0625,.0025)

fig=MapContourf(lon_targcyc,lat_targcyc,sigma_avg_calc_day,
                levels=clevs,
                figsize=(14,5),
                extend='both',
                title="Eady inverse time scale calc $(days^{-1}) \sim {|f| \over \sqrt{R_i}}$")
maxf=np.maximum(np.sqrt(c_rossby_nearest*2.*btp_nearest),absf)
sigma_maxf_calc_day=(maxf/sqrtgrate_nearest)*86400.
clevs=arange(.005,.0625,.0025)

#fig=MapContourf(lon_targcyc,lat_targcyc,sigma_maxf_calc_day,
#                levels=clevs,
#                figsize=(14,5),
#                extend='both',
#                title="Eady inverse time scale calc maxf $(days^{-1}) \sim {|f| \over \sqrt{R_i}}$")

-1e+30 -1e+30

In [43]:

tulloch_eady=Image(filename=basedir+'/steering/tulloch_eady.png')
display(tulloch_eady)

Alternate Sigma Calculation

$\sigma_{vi} = {f \over \sqrt{R_i}}\quad = \quad {f \over \sqrt{{f^2N^2 \over m^4}}}\quad = \quad {{\cancel f m^2} \over \cancel f N}$

In [44]:

clevs=arange(.005,.0625,.0025)

fig=MapContourf(lon_targcyc,lat_targcyc,sigma_avg_calc_day,
                levels=clevs,
                figsize=(14,5),
                extend='both',
                title="Old Eady inverse time scale calc $(days^{-1}) \sim {|f| \over \sqrt{R_i}}$")

ncalt = Dataset(basedir+'/steering/g.e11.GIAF.T62_gx1v6.steer.007.pop.h.nmonth1.0249-02.nc')
#print(nc.variables)
sigma_avg1 = ncalt.variables['SIGMA_AVG'][0,:,:]
sigma_avg1_nearest = pyresample.kd_tree.resample_nearest(orig_def, sigma_avg1, \
        targ_def, radius_of_influence=500000, fill_value=None)
sigma_avg1_monthavg_day=sigma_avg1_nearest*86400.
clevs=arange(.005,.0625,.0025)

fig=MapContourf(lon_targcyc,lat_targcyc,sigma_avg1_monthavg_day,
                levels=clevs,
                figsize=(14,5),
                extend='both',
                title="New Eady inverse time scale calc using alt sigma $(days^{-1}) \sim  {{m^2} \over N}$")

In [45]:

bates_len_vel=Image(filename=basedir+'/steering/bates_len_vel.jpg')
display(bates_len_vel)

Zonal Eddy phase speed:

$\cancel{c = - \beta * {L_r^2}}$ $L_r$ too high at equator

$c = - \beta * L_{eddy}^2$

In [46]:

#beta
clevs=arange(1.e-14,25.e-14,1e-14)
norm=matplotlib.colors.Normalize(),
fig12=MapContourf(lon_targcyc,lat_targcyc,btp_nearest,
                    levels=clevs,
#                    norm=norm,
                    addzonal=True,
                    figsize=(14,5),
                    extend='both',
                    title="Beta (1/cm*s)")


clevs=np.logspace(1,15,30)
norm=LogNorm()
l_eddy_calc_sq=l_eddy_calc*l_eddy_calc
fig8=MapContourf(lon_targcyc,lat_targcyc,l_eddy_calc_sq,
                    addzonal=True,
                    levels=clevs,
                    norm=norm,
                    figsize=(14,5),
                    extend='both',
                    title="$L_{eddy}^2 (cm^2)$ log scale")

c=btp_nearest*l_eddy_calc*l_eddy_calc # units=1/(cm s) * cm^2 = cm/s

#clevs=np.linspace(-10,30,10)
clevs=arange(-10,30,1)
norm=LogNorm()
fig12=MapContourf(lon_targcyc,lat_targcyc,c,
                    addzonal=True,
                    levels=clevs,
#                    norm=norm,
                    figsize=(14,5),
                    contourlines=True,
                    extend='both',
                    title="eddy phase speed c = beta*leddy**2 (cm/s) linear scale(NOTE I'm plotting positive c here - changed to negative after this plot")

c=-1.*c

Hughes Phase Speed (cm/s) from Tulloch Marshall Smith '09

In [47]:

hughes_c=Image(filename=basedir+'/steering/Hughes_phase_speed.png')
display(hughes_c)

Eddy Velocity:

$u_{rms} = alpha*\sigma*L_{eddy}$

In [48]:

############################
##### PLOT L_Eddy
############################

cm2m=.01
alpha=5.

sigma_calc=(absf/sqrtgrate_nearest)
#sigma_calc=(maxf/sqrtgrate_nearest)

u_rms_leddy = alpha*sigma_calc*l_eddy_calc    # units = 1*1/s*cm = cm/s
u_rms_leddy_mps=u_rms_leddy*cm2m

clevs=np.logspace(0,6,20)
cm2m=.01
l_eddy_calc_m=l_eddy_calc*cm2m
norm=LogNorm()
fig9=MapContourf(lon_targcyc,lat_targcyc,l_eddy_calc_m,
                    addzonal=True,
                    levels=clevs,
                    norm=norm,
                    figsize=(14,5),
#                    extend='both',
                    title="$L_{eddy} (m)$ log scale")
pylab.xlim([0,3e5])
#################  PLOT u_rms with Leddy #######################
u_rms_leddy_mps=u_rms_leddy*.01
clevs=np.linspace(0,.4,40)
fig99=MapContourf(lon_targcyc,lat_targcyc,u_rms_leddy_mps,
                    addzonal=True,
                    levels=clevs,
                    figsize=(14,5),
                    extend='both',
                    title="$u_{rms} ({m \over s})$ (alpha=%i) %i * old $\sigma*l_{eddy}$"%(alpha,alpha))
pylab.xlim([0,.2])

sigma_calc_new=sigma_avg1_nearest
u_rms_leddy_new = alpha*sigma_calc_new*l_eddy_calc    # units = 1*1/s*cm = cm/s
u_rms_leddy_new_mps=u_rms_leddy_new*cm2m

clevs=np.linspace(0,.4,40)
fig99=MapContourf(lon_targcyc,lat_targcyc,u_rms_leddy_new_mps,
                    addzonal=True,
                    levels=clevs,
                    figsize=(14,5),
                    extend='both',
                    title="$u_{rms} ({m \over s})$ (alpha=%i) %i * alt $\sigma*l_{eddy}$"%(alpha,alpha))
pylab.xlim([0,.2])

Out[48]:

(0, 0.2)

In [49]:

bates_urms=Image(filename=basedir+'/steering/bates2013-urms.png')
display(bates_urms)

In [50]:

#################  PLOT u_rms^2  convert m^2/s^2 #######################
u_rms_leddy_mps_sq=u_rms_leddy_mps*u_rms_leddy_mps
#clevs=np.linspace(0,.06,50)
clevs=arange(0,.06,.001)
#clevs=np.logspace(-5,1,30)
norm=LogNorm()
fig99=MapContourf(lon_targcyc,lat_targcyc,u_rms_leddy_mps_sq,
                    addzonal=True,
                    levels=clevs,
#                    norm=norm,
                    figsize=(14,5),
                    extend='both',
                    title="$u_{rms}^2 ({m^2 \over s^2})$ (alpha=%i)   %i * old $\sigma*l_{eddy}$ "%(alpha,alpha))
pylab.xlim([0,.06])

u_rms_leddy_new_mps_sq=u_rms_leddy_new_mps*u_rms_leddy_new_mps
#clevs=np.linspace(0,.06,50)
clevs=arange(0,.06,.001)
#clevs=np.logspace(-5,1,30)
norm=LogNorm()
fig99=MapContourf(lon_targcyc,lat_targcyc,u_rms_leddy_new_mps_sq,
                    addzonal=True,
                    levels=clevs,
#                    norm=norm,
                    figsize=(14,5),
                    extend='both',
                    title="$u_{rms}^2 ({m^2 \over s^2})$ (alpha=%i)   %i * Alt $\sigma*l_{eddy}$ "%(alpha,alpha))
pylab.xlim([0,.06])

Out[50]:

(0, 0.06)

In [51]:

urms_bates=Image(filename=basedir+'/steering/Bates_urms_sq.jpg')
display(urms_bates)

In [52]:

#################  PLOT (u-c)^2  convert m^2/s^2 #######################
umean10_nearest = pyresample.kd_tree.resample_nearest(orig_def, umean10,
                targ_def, radius_of_influence=500000, fill_value=None)

umean10_mps=umean10_nearest*1.e-2

clevs=np.linspace(-60,60,40)
fig11=MapContourf(lon_targcyc,lat_targcyc,umean10_nearest,
                    addzonal=True,
                    figsize=(14,5),
                    levels=clevs,
                    extend='both',
                    title="U zonal velocity (cm/s)")
pylab.xlim([-20.,20.])

# get 2D versions of the lat and lon variables add longitude start point here!
lon_ecco_orig, lat_ecco_orig = np.meshgrid(lon_ecco[:], lat_ecco[:])

orig_ecco_def = pyresample.geometry.GridDefinition(lons=lon_ecco_orig, lats=lat_ecco_orig)

ecco_u_nearest = pyresample.kd_tree.resample_nearest(orig_ecco_def, ecco_u, \
        targ_def, radius_of_influence=500000, fill_value=None)
ecco_u_nearest_masked=np.ma.masked_where(ecco_u_nearest<-10, ecco_u_nearest, copy=True)

#print ecco_u_cyc_masked[:,0]

clevs=np.linspace(-60,60,40)
#ecco_u_nearest = pyresample.kd_tree.resample_nearest(orig_def, ecco_u, \
#        targ_def, radius_of_influence=500000, fill_value=None)
fig11=MapContourf(lon_targcyc,lat_targcyc,ecco_u_nearest_masked*100,
                    addzonal=True,
                    figsize=(14,5),
                    levels=clevs,
#                    extend='both',
                    title="ECCO U zonal velocity (cm/s)")
pylab.xlim([-20.,20.])

Out[52]:

(-20.0, 20.0)

In [53]:

clevs=-np.logspace(3,-3,20)
print clevs
norm=matplotlib.colors.SymLogNorm(linthresh=1e-4, vmin=-1000, vmax=-.001)
fig12=MapContourf(lon_targcyc,lat_targcyc,c,
                    addzonal=True,
                    levels=clevs,
                    norm=norm,
                    figsize=(14,5),
#                    extend='both',
                    title="eddy phase speed c = beta*leddy**2 (cm/s) log scale")
pylab.xlim([-200.,0.])

uminusc=umean10_nearest-c  #units cm/s
#clevs=np.logspace(-2,3,10)
clevs=np.linspace(0,25,25)
#norm=LogNorm()
fig13=MapContourf(lon_targcyc,lat_targcyc,uminusc,
                    addzonal=True,
#                    norm=norm,
                    levels=clevs,
                    figsize=(14,5),
#                    extend='both',
                    title="U-c (cm/s) (max = %i/min = %i) (c calculated using Leddy)"%(uminusc.max(),uminusc.min()))
pylab.xlim([0,25])

[ -1.00000000e+03  -4.83293024e+02  -2.33572147e+02  -1.12883789e+02
  -5.45559478e+01  -2.63665090e+01  -1.27427499e+01  -6.15848211e+00
  -2.97635144e+00  -1.43844989e+00  -6.95192796e-01  -3.35981829e-01
  -1.62377674e-01  -7.84759970e-02  -3.79269019e-02  -1.83298071e-02
  -8.85866790e-03  -4.28133240e-03  -2.06913808e-03  -1.00000000e-03]

Out[53]:

(0, 25)

In [54]:

bates_uminusc=Image(filename=basedir+'/steering/bates2013_uminusc.png')
display(bates_uminusc)

In [55]:

uminusc_mps_sq=uminusc*uminusc*1.e-4   #units cm/s*cm/s*m/100cm*m/100cm
clevs=np.linspace(0,.06,41)
fig14=MapContourf(lon_targcyc,lat_targcyc,uminusc_mps_sq,
                    addzonal=True,
                    levels=clevs,
                    figsize=(14,5),
#                    extend='both',
                    title=" ${(U-c)}^2 {(m^2/s^2)}$ max = %i/min = %i (c calculated using Leddy) "%(uminusc_mps_sq.max(),uminusc_mps_sq.min()))

In [56]:

bates_uminuscsq=Image(filename=basedir+'/steering/Bates_u_minus_c_sq.jpg')
display(bates_uminuscsq)

In [57]:

uminusc_sq_ovr_urms_sq = uminusc_mps_sq/u_rms_leddy_mps_sq             #units cm/s /  cm/s  
uminusc_sq_ovr_urms_sq_masked=np.ma.masked_where(uminusc_sq_ovr_urms_sq >30., uminusc_sq_ovr_urms_sq, copy=True)
clevs=arange(0,5,.1)
#clevs=np.logspace(-3,3,21)
norm=LogNorm()
fig14=MapContourf(lon_targcyc,lat_targcyc,uminusc_sq_ovr_urms_sq,                
                    addzonal=True,
                    levels=clevs,
#                    norm=norm,
                    figsize=(14,5),
                    extend='both',
                    title="${(U-c)}^2 \over u_{rms}^2$")
pylab.xlim([0,8])

uminusc_sq_ovr_urms_new_sq = uminusc_mps_sq/u_rms_leddy_new_mps_sq             #units cm/s /  cm/s  
clevs=arange(0,5,.1)
#clevs=np.logspace(-3,3,21)
norm=LogNorm()
fig14=MapContourf(lon_targcyc,lat_targcyc,uminusc_sq_ovr_urms_new_sq,                
                    addzonal=True,
                    levels=clevs,
#                    norm=norm,
                    figsize=(14,5),
                    extend='both',
                    title="${(U-c)}^2 \over u_{rms}^2$ using Alt sigma")
pylab.xlim([0,8])


clevs=np.logspace(-3,4,50)
norm=LogNorm()
fig14=MapContourf(lon_targcyc,lat_targcyc,uminusc_sq_ovr_urms_new_sq,                
                    addzonal=True,
                    levels=clevs,
                    norm=norm,
                    figsize=(14,5),
                    extend='both',
                    title="${(U-c)}^2 \over u_{rms}^2$ max=%i/min=%i (using Alt sigma) log scale"%(uminusc_sq_ovr_urms_new_sq.max(),uminusc_sq_ovr_urms_new_sq.min()))
pylab.xlim([0,100])

Out[57]:

(0, 100)

In [58]:

bates_uminusc_sq_ovr_urms_sq =Image(filename=basedir+'/steering/Bates_suppress.jpg')
display(bates_uminusc_sq_ovr_urms_sq )

$K=u_{rms}∗L_{mix}$

$L_{mix} = \Gamma * L_{eddy} * Suppression$

$Suppression= {1 \over (1 + b1 * |u_{mean} - c|^2 /u_{rms (z=0)}^2 )}$

$L_{mix} = {\Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - c|^2 /u_{rms}^2 (z=0)}$

In [59]:

cm2m=.01
u_rms_leddy_mps=u_rms_leddy*cm2m
u_rms_leddy_mps_sq=u_rms_leddy_mps*u_rms_leddy_mps
gamma=.35
b1=4.0
supp=1./(1+b1*uminusc_sq_ovr_urms_sq)
clevs=np.linspace(0,1,30)
clevs=arange(0,1,.05)
norm=matplotlib.colors.Normalize()
fig=MapContourf(lon_targcyc,lat_targcyc,supp,
                    addzonal=True,
                    norm=norm,
                    levels=clevs,
                    figsize=(14,5),
                    title="Suppression factor ${1 \over (1 + b1 * |u_{mean} - c|^2 /u_{rms (z=0)}^2)}$ using old sigma")

pylab.xlim([0,.8])
supp_new=1./(1+b1*uminusc_sq_ovr_urms_new_sq)
clevs=np.linspace(0,1,20)
norm=matplotlib.colors.Normalize()
fig=MapContourf(lon_targcyc,lat_targcyc,supp_new,
                    addzonal=True,
                    norm=norm,
                    levels=clevs,
                    figsize=(14,5),
                    title="Suppression factor ${1 \over (1 + b1 * |u_{mean} - c|^2 /u_{rms (z=0)}^2)}$ using Alt sigma")
pylab.xlim([0,.8])

Out[59]:

(0, 0.8)

In [60]:

bates_suppression =Image(filename=basedir+'/steering/bates2013_suppression.png')
display(bates_suppression )

In [61]:

gammax5=gamma*5
lmix_new_m=gammax5*l_eddy_calc_m*supp_new
clevs=np.logspace(0,6,20)
norm=LogNorm()
fig9=MapContourf(lon_targcyc,lat_targcyc,l_eddy_calc_m,
                    addzonal=True,
                    levels=clevs,
                    norm=norm,
                    figsize=(14,5),
#                    extend='both',
                    title="$L_{eddy} (m)$ log scale")
pylab.xlim([0,3e5])



clevs=np.linspace(0,100e3,10)
norm=matplotlib.colors.Normalize()
fig=MapContourf(lon_targcyc,lat_targcyc,lmix_new_m,
                    addzonal=True,
                    norm=norm,
                    levels=clevs,
                    figsize=(14,5),
                    title="lmix (m) = gamma*leddy*suppression (using Alt sigma ... also used gamma x 5  (5*.35))")
#pylab.xlim([0,10e3])

clevs=np.linspace(0,.4,40)
fig99=MapContourf(lon_targcyc,lat_targcyc,u_rms_leddy_new_mps,
                    addzonal=True,
                    levels=clevs,
                    figsize=(14,5),
                    extend='both',
                    title="$u_{rms} ({m \over s})$ (alpha=%i) %i * alt $\sigma*l_{eddy}$"%(alpha,alpha))
pylab.xlim([0,.2])

k_new=u_rms_leddy_new_mps*lmix_new_m
clevs=np.linspace(0,10e3,40)
norm=matplotlib.colors.Normalize()
fig=MapContourf(lon_targcyc,lat_targcyc,k_new,
                    addzonal=True,
                    norm=norm,
                    levels=clevs,
                    figsize=(14,5),
                    title="k=urms*lmix (using Alt sigma for calculation of urms and gamma * 5 for calculation of lmix")
pylab.xlim([0,10e3])

Out[61]:

(0, 10000.0)

In [62]:

urms_bates2=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-002.jpg')
urms_bates3=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-003.jpg')
urms_bates4=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-004.jpg')
urms_bates5=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-005.jpg')
urms_bates6=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-006.jpg')
urms_bates7=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-007.jpg')
urms_bates8=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-008.jpg')
urms_bates9=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-009.jpg')
urms_bates10=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-010.jpg')
urms_bates11=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-011.jpg')
display(urms_bates3,urms_bates4,urms_bates5,urms_bates6,urms_bates7,urms_bates8,urms_bates9,urms_bates10,urms_bates11)

CESM $U_{mean} = U_{resolved}$ integraged over {10,50,100,200,500} meters

In [63]:

umean_cesm=Image(filename=basedir+'/steering/umean.jpg')
display(umean_cesm)

In [64]:

sigma_avg_calc=(absf/sqrtgrate_nearest)
rhines_km=sigma_avg_calc/btp_nearest*cm2km
print rhines_km.max(),rhines_km.min()
clevs=[10,20,30,40,50,60,80,100,150,225]
clevs=arange(0,225,10)
fig=MapContourf(lon_targcyc,lat_targcyc,rhines_km,
                addzonal=True,
                levels=clevs,
                figsize=(14,5),
                extend='both',
                title="Rhines Scale $(km) \sim {Sigma \over B}$")

-- --

In [65]:

    popdata = Dataset(tfile, 'r')
    print tfile
    print popdata
    shear_max = popdata.variables('SHEAR_MAX')

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-65-bc632d084a56> in <module>()
----> 1 popdata = Dataset(tfile, 'r')
      2 print tfile
      3 print popdata
      4 shear_max = popdata.variables('SHEAR_MAX')

NameError: name 'tfile' is not defined

In []:

from ctd import plot_section
from oceans.colormaps import cm
#from oceans.datasets import etopo_subset
def laplace_X(F, M):
    """1D Laplace Filter in X-direction.
    http://www.trondkristiansen.com/wp-content/uploads/downloads/
    2010/09/laplaceFilter.py"""

    jmax, imax = F.shape

    # Add strips of land.
    F2 = np.zeros((jmax, imax + 2), dtype=F.dtype)
    F2[:, 1:-1] = F
    M2 = np.zeros((jmax, imax + 2), dtype=M.dtype)
    M2[:, 1:-1] = M

    MS = M2[:, 2:] + M2[:, :-2]
    FS = F2[:, 2:] * M2[:, 2:] + F2[:, :-2] * M2[:, :-2]

    return np.where(M > 0.5, (1 - 0.25 * MS) * F + 0.25 * FS, F)


def laplace_Y(F, M):
    """1D Laplace Filter in Y-direction.
    http://www.trondkristiansen.com/wp-content/uploads/downloads/
    2010/09/laplaceFilter.py"""

    jmax, imax = F.shape

    # Add strips of land.
    F2 = np.zeros((jmax + 2, imax), dtype=F.dtype)
    F2[1:-1, :] = F
    M2 = np.zeros((jmax + 2, imax), dtype=M.dtype)
    M2[1:-1, :] = M

    MS = M2[2:, :] + M2[:-2, :]
    FS = F2[2:, :] * M2[2:, :] + F2[:-2, :] * M2[:-2, :]

    return np.where(M > 0.5, (1 - 0.25 * MS) * F + 0.25 * FS, F)


def laplace_filter(F, M=None):
    """Laplace filter a 2D field with mask.  The mask may cause laplace_X and
    laplace_Y to not commute.  Take average of both directions.
    http://www.trondkristiansen.com/wp-content/uploads/downloads/
    2010/09/laplaceFilter.py"""

    if not M:
        M = np.ones_like(F)

    return 0.5 * (laplace_X(laplace_Y(F, M), M) +
                  laplace_Y(laplace_X(F, M), M))


def pop_subset(var='',llcrnrlon=None, urcrnrlon=None, llcrnrlat=None, urcrnrlat=None, tfile='', smoo=False, subsample=False):

### Read Pop file
    popdata = Dataset(tfile, 'r')
    popvar = popdata.variables[var][0,:,:,:]
    
    print popvar.shape
    popvar=np.rollaxis(popvar, 1)
    print popvar.shape
    popvar=np.rollaxis(popvar, 2,1)
    print popvar.shape

### Get lats lons for original and target grids for remapping 
### Process input lons for pyregrid (-180:180) and output lons for matplotlib idiosyncracies

    lon_orig = popdata.variables["TLONG"][:]
    lat_orig = popdata.variables["TLAT"][:]
    lev_targ = popdata.variables["z_t"][:]
# setup for using pyresample to interpolate plot
# make sure lon is in range of -180:180 for pyregrid
    lon_orig = (lon_orig + 180) % (360) - 180 

# open dataset for output grid get latout,lonout
    fc = Dataset(basedir+'/steering/co2flux_fossil_1751-2006-monthly_0.9x1.25_c20100204.nc')
    lon= fc.variables['lon'][:]
    lat = fc.variables['lat'][:]
    fc.close()

# get 2D versions of the lat and lon variables add longitude start point here!
    lon_targ, lat_targ = np.meshgrid(lon, -1.*lat)

# add cyclic point to output lat,lon arrays for plot routines
    lon_targcyc = np.concatenate((lon_targ, lon_targ[:, 0:1]), axis=1)
    lat_targcyc = np.concatenate((lat_targ, lat_targ[:, 0:1]), axis=1)

# make sure lon is in range of -180:180 for pyregrid
    lon_targcyc = (lon_targcyc + 180) % (360) - 180

#make longitudes monotonic for contouring add longitude period to end
    lon_targ = np.where(np.greater_equal(lon_targcyc,min(lon_targcyc[:,0])),lon_targcyc-360,lon_targcyc)+360.
    endlon=lon_targ.shape[1]-1
    lon_targ[:,endlon]=lon_targ[:,endlon]+360.
    
    lat_targ=lat_targcyc

#When all the necessary inputs are read in memory, that's where we start using pyresample's magic.
#Create a pyresample object holding the origin (POP) grid:

    orig_def = pyresample.geometry.GridDefinition(lons=lon_orig, lats=lat_orig)

#Create another pyresample object for the target (curvilinear) grid:
    targ_def = pyresample.geometry.GridDefinition(lons=lon_targcyc, lats=lat_targcyc)
    print lon_orig.shape,lat_orig.shape,lon_targ.shape,lat_targ.shape,popvar.shape
    popvar_remap = pyresample.kd_tree.resample_nearest(orig_def, popvar, \
        targ_def, radius_of_influence=500000, fill_value=None)


    # Select data subset.
maskx = np.logical_and(lon_targ >= llcrnrlon, lon_targ <= urcrnrlon)
masky = np.logical_and(lat_targ >= llcrnrlat, lat_targ <= urcrnrlat)
maskz = np.array([z in depths for z in depth])

    
    ### set region limits

###res = get_indices(llcrnrlat, urcrnrlat, llcrnrlon, urcrnrlon, lon_targ, lat_targ)

###lon, lat = np.meshgrid(lons[res[0]:res[1]], lats[res[2]:res[3]])

####popvar_clip = ncdataset.variables[var][int(res[2]):int(res[3]),int(res[0]):int(res[1])]

#    if smoo:
#        popvar_remap = laplace_filter(popvar_remap, M=None)

    if subsample:
        lon_targ, lat_targ, popvar_remap = lon_targ[::subsample], lat_targ[::subsample], popvar_remap[::subsample]

    return lon_targ, lat_targ, lev_targ, popvar_remap

In []:

from __future__ import division

import numpy as np
import matplotlib.pyplot as plt

from pandas import Panel4D, Panel
from netCDF4 import Dataset, num2date
from ..ff_tools import get_profile, wrap_lon180, wrap_lon360

#
def map_limits(m):
    lons, lats = wrap_lon360(m.boundarylons), m.boundarylats
    boundary = dict(llcrnrlon=min(lons),
                    urcrnrlon=max(lons),
                    llcrnrlat=min(lats),
                    urcrnrlat=max(lats))
    return boundary

longitude = -25.5
llcrnrlon=longitude
urcrnrlon=longitude

llcrnrlon, urcrnrlon = map(wrap_lon360, (llcrnrlon, urcrnrlon))

lon = wrap_lon360(nc.variables.pop('lon')[:])
lat = nc.variables.pop('lat')[:]
depth = nc.variables.pop('depth')[:]
times = nc.variables.pop('time')
times = num2date(times[:], times.units, calendar='365_day')
times = [time.strftime('%b') for time in times]

longitude = -25.5
llcrnlon=longitude

tfile=basedir+'/steering/g.e11.GIAF.T62_gx1v6.steer.007.pop.h.nday1.0249-01-01.nc'
nc= Dataset(tfile, 'r')
lon = wrap_lon360(nc.variables.pop('lon')[:])
print lon
    # Select data subset.
#maskx = np.logical_and(lon >= llcrnrlon, lon <= urcrnrlon)
#masky = np.logical_and(lat >= llcrnrlat, lat <= urcrnrlat)
#maskz = np.array([z in depths for z in depth])


levels=slice(0,20)
depths = [0, 10, 20, 30, 50, 75, 100, 125, 150, 200, 250, 300, 400, 500,
              600, 700, 800, 900, 1000, 1100, 1200, 1300, 1400, 1500, 1750,
              2000, 2500, 3000, 3500, 4000, 4500, 5000, 5500, 6000, 6500, 7000,
              7500, 8000, 8500, 9000][levels]
print depths
var='temperature'
v = dict(temperature='t', dissolved_oxygen='o', salinity='s',
             oxygen_saturation='O', apparent_oxygen_utilization='A',
             phosphate='p', silicate='p', nitrate='n')
start = '%s_' % v[var]
print start
variables = dict()
for variable in nc.variables.keys():
    if variable.startswith(start):
        subset = nc.variables[variable][..., maskz, masky, maskx]
        data = Panel4D(subset, major_axis=lat, minor_axis=lon,
                        labels=np.atleast_1d(times),
                        items=np.atleast_1d(depth))
        variables.update({d[variable]: data})
print variables

In []:

#dataset = dataset.clip(27.5, 37.5)
#levels = np.arange(dataset.min().min(), dataset.max().max(), 0.02)
#dataset.lat = dataset.columns.values.astype(float)
#dataset.lon = len(dataset.lat) * [longitude]

tfile=basedir+'/steering/g.e11.GIAF.T62_gx1v6.steer.007.pop.h.nday1.0249-01-01.nc'
longitude=-25.5
lons, lats, levs, shear = pop_subset(var='SHEAR_MAX',
                               llcrnrlon=longitude,
                               urcrnrlon=longitude,
                               llcrnrlat=-70,
                               urcrnrlat=70,
                               smoo=True,
                               tfile=tfile)
print shear.shape

print shear[:,240,:].shape
print levs.shape
fig, ax = plt.subplots()
latm, levm = np.meshgrid(lats[:,240], levs[::-1])
cs = ax.pcolormesh(latm, levm, shear[:,240,:].transpose())

fig, ax, cbar = plot_section(shear, cmap=ocm.avhrr, marker=None,
                             levels=levels, figsize=(10, 4))
X = ax.get_xticks()
offset = 1.01
ax.set_xlim(X.min() * offset, X.max() * offset)
new_labels = np.linspace(dataset.lat.min(), dataset.lat.max(), len(X))
new_labels = [u'%i\u00B0' % num for num in new_labels]
ax.set_xticklabels(new_labels)
ax.set_xlabel(u'WOA09 Salinity Section at Longitude %3.1f\u00B0' % longitude)
ax.set_ylabel('Depth [m]')

axin = inset_axes(ax, width="40%", height="40%", loc=4)
inmap = Basemap(projection='ortho', lon_0=longitude, lat_0=0,
                ax=axin, anchor='NE')
inmap.bluemarble()
_ = inmap.plot(dataset.lon, dataset.lat, 'r', latlon=True)

POP Steering parameterization

Eddy diffusivity implementation by Bates, Tulloch, Marshall, and Ferrari

POP Steering Parameterization

CAM Implementation of Mesoscale Eddy Diffusivity in Terms of Mixing Length Theory

MICHAEL BATES, ROSS TULLOCH, JOHN MARSHALL, AND RAFFAELE FERRARI

Steering Level status

Todo

Update the matplotlib configuration parameters:

mpl.rcParams.update({'font.size': 18, 'font.family': 'serif'})

Implementation Notes on Diffusivity with Steering Level Suppression

Eddy Length Scales:

Eddy Velocity:

An alternate derivation of $\sigma_{vi}$ is:

Implementation notes

For Level K

CESM First Baroclinic Wave Speed $c_r \cancel {= {1\over\pi} \int_{-2000}^{-100} N(z)dz}$ Now it is over full model depth

First Baroclinic Rossby Radius Chelton vs CESM

CESM Rossby deformation radius = $L_{eddy}$

$\sigma$ Eady growth rate

Alternate Sigma Calculation

Zonal Eddy phase speed:

Hughes Phase Speed (cm/s) from Tulloch Marshall Smith '09

Eddy Velocity:

$K=u_{rms}∗L_{mix}$

CESM $U_{mean} = U_{resolved}$ integraged over {10,50,100,200,500} meters